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Improvements of a congruential method for finding the exact solutions of systems of linear equa- 
tions with rational coefficients are described. Typical execution times on the CDC 1604-A are given, 
as well as the Fortran program. 
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1 . Introduction 

An algorithm for computing the exact solutions of linear equations with rational coefficients, 
and its computer implementation, were described in [l]. 1 The basic idea of the algorithm is to 
convert the original system of equations into a system of congruences modulo various primes 
Pi, and combining the solutions by a procedure suggested by the Chinese Remainder Theorem. 
This process is continued until the sequence of solutions modulo p\p>> . . . Pk, fc=l, 2, . . . 
converges to the true solution. The major part of the computation is performed in single precision. 
See also Newman [4] who used the method for computing the inverses of ill-conditioned matrices, 
and Knuth [3, p. 256], who remarked that for ill-conditioned matrices the procedure "gives a 
method for obtaining the true answers in less time than any known method can produce reliable 
approximate answers!" 

Our program provides for a final substitution check for verifying convergence, i.e., verifying 
that the computed values satisfy the original system. It is applied after two successive iterations 
produce no change in the solution vectors. If so implemented, and assuming sufficient memory 
space, the method produces the exact solution space for any solvable linear system with rational 
coefficients. Moreover, the algorithm is effective, in the sense that the exact solution space is 
produced within a reasonable time for systems that are not too large. It cannot end up with a 
wrong or no solution, as additional iterations are made should the substitution check fail. 

The implementation of the algorithm has since been improved. As a result, the computation 
of the 9 independent exact solutions of a system of 111 homogeneous equations in 120 unknowns 
of rank 111 with integral coefficients in the range [—2180, 2568] which took 60 min on a CDC 
1604^-A by the old method, now takes only 19 min on the same computer. The program is still 
in the form of standard Fortran subroutines, and no pains were taken to write a particularly eco- 
nomic program. 

In the sequel it is assumed that the reader is reasonably familiar with the essential features 
of [1]. 
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2. The Improvements 

There are four main improvements: Use of larger primes; using the Cantor representation 
for constructing the solution mod p\pz . . . pu after the kth iteration (by the Chinese Remainder 
Theorem); a logical simplification in determining whether to retain or discard a new prime; and 
using a faster algorithm for finding the greatest common divisor (g.c.d.) of the components of 
any solution vector. These improvements will now be described one for one. 

(i) Use of larger primes. In [1], primes of the order of 10 7 , just less than half the machine 
word length, were employed. Primes of the order of 10 14 , just less than one machine word are now 
used. Such primes were supplied to us through the courtesy of Herschel F. Smith from IBM. 
This saves about half the number of iterations (primes). In the above example, it reduced the 
number of iterations from six to four. Of course, we could have used the squares of primes of the 
order of 10 7 to obtain the same effect. However, the probability that a prime divides any of the 
principal minors of a matrix decreases with the size of the prime. Hence it is of some advantage 
to use the largest primes just less than a machine word. 

(ii) Use of the Cantor representation. Let mi, . . ., m s be odd positive integers, relatively 
prime in pairs. Any number N in the range 



m\ni2 . . . m, s 



TfllHfe • • • ™>s 

<N^ 



is uniquely representable by the Cantor representation, also referred to, by several authors, as 
the mixed radix representation: 

N=a,\ + a%m\ -\- a^mxirh + . . . H- a s m\ /n 2 . . . m s -\, 
where 



m, 

Y 



< a% 



Y 



i = l, 2, ... ,5. 



The Cantor representation of a number given by its residues mod rrtu i— 1, 2, . . ., s, can 
be determined by computations in which all numbers occurring have absolute value not exceeding 
max mi, in the present case single precision. See e.g., [2], The solution check of [1] is now per- 
formed when for any fixed k, a^— for allcomponents of the solution vectors. 

Thus, the main computation of the solution vectors is now performed in single precision, 
resulting in a large time saving. Multiple precision is required only in the final conversion of the 
Cantor representation to decimal representation, and in the — optional— subsequent computation 
of the g.c.d. of the components of each solution vector. 

(hi) A logical simplification. On pp. 110-111 of [1], a method was described which guaran- 
tees that each triangularization converges to the same largest nonsingular submatrix of the coeffi- 
cient matrix A. For this, a lexicographic ordering of rows and columns of A had to be checked. This 
has now been changed as follows: The triangularization of A (mod p\) induces a certain interchange 
of rows and columns in A. If p is the rank of A mod p u let i u . . .,i p ,/i, . . ., j P be the rows and 
columns of A appearing in the triangularized principal submatrix of order p. The rows and columns 
of A itself are now reordered, so that i u . . ., i p , j\, . . . , j p appear as its first p rows and col- 
umns. Call this permuted matrix B. In each subsequent iteration mod pi (i> 1), the triangulariza- 
tion is performed without changing any of the first p rows and columns in B. In the triangularized 
matrix, denote by pi the order of the largest nonsingular principal submatrix. If p* < p, pt must be 
discarded; if pi > p, the primes p\ , . . . , pt-i must be discarded and a new matrix B with a greater 
rank is formed from A; if pi = p, which is the normal case, the solution mod pi is used to determine 
the coefficients a\ of the Cantor representation of all components of the solution vectors. 

This modification does not save much time. In practice it only saves logical operations and 
manipulations, as it almost never happens that a given large prime divides a nonzero principal 
minor. However, the logical simplification results in a more compact and elegant program. 
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(iv) Computation of the g.c.d. The last step in [1] was to find the g.c.d. of the components of 
each of the solution basis vectors and to divide them out, so as to obtain primitive solutions, i.e., 
solutions such that the g.c.d. of the components of each solution vector is unity. Instead of com- 
puting the g.c.d. by the Euclidean Algorithm, i.e., by a series of multiple precision divisions, it 
is found by a series of shifts, using an economic algorithm of Stein [5], which, according to Knuth 
[3, p. 297], was previously given by Silver and Terzian: 

Let t be any positive integer stored in binary form in a computer register. The highest power 
of 2 dividing t can easily be determined by shifting t to the right until its least significant nonzero 
bit is located in the least significant position of the register. This results in an odd integer t' = 2~ k t. 
We shall designate this operation by shift, i.e., t' = shifty. 

Let a, b be two positive integers, and let 

a = shift {a) = 2~ k a, b = shift (b) = 2~ l b. 

Let m=min (k, I). The algorithm now proceeds as follows: 

«i = shift |a — 6o|, b\ = min (a , bo) 
a 2 = shift \tii — 6i|, 6 2 — min (a u bi) 



a n — shift |a„_i— 6„_i|, b n = min (a n -i, b n -i) 
rtn+i = shift \a n — b n \. 

The procedure terminates when a n = b n , and then (a, b) = 2 m a n . 

In our example, a multiple precision g.c.d. was obtained for each of the 9 solution vectors. 
The computation of these g.c.d. and dividing them out now takes 2.5 min instead of the previous 
7.5 min. 

Finally, we should remark that the substitution check, which is also a multiple-precision 
operation, takes 6 min for our example. It is extremely rare that a substitution check fails when 
large primes are used, and therefore it seems reasonable to dispose of it in general. In fact, no 
substitution check is normally made in any of the conventional iterative schemes for solving 
linear equations — including those for which convergence is not guaranteed a priori even when the 
system is known to be solvable. 

On the other hand, it is easy to fabricate a failing case. For simplicity assume that there is 
only one solution vector V = (v\, v 2 , . • ., vt) , whose Cantor representation is 

Vi=an+ CL12P1 + a&piP2+ . . . + ciikPiP2 . . . Pk-i + • • ■ + auPiP2 • • • Ps-u 1 ^ i ^ t, 

where Pi, /?2, • - • are distinct primes. For obtaining a failing case, we simply choose a^ = 0, 
l^i^t, and ay, ;V k arbitrarily, with the only condition that an ^ for some l> k and for 
some i. 

A vector V oi this form can clearly also be characterized as follows. Suppose that V k ~ x = (v k ~\ 
v k ~ l , . . ., v k ~ x ) is the solution mod pip> . . . Pk-u i.e., 



Then 

if and only if 



Vi = v k ~ l (mod pip 2 . . . P/c-i), 1 ^ i ^ t. 

cLi k = , 1 ^ i ^ t, an t^ 0, / > k for some i 

Vi = z;f _1 (mod p^) , 1 ^ i ^ t, Vi ¥" v^~ l for some i. 



Disposing of the substitution check and the g.c.d. computation, solution time in our example 
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reduces to 10 min. By comparison, the approximate computation of a well-conditioned 111 X 111 
system by standard numerical floating number techniques takes about 3 min on the CDC 1604 -A. 

3. Machine Program 

The program given below is the full version, i.e., it includes computation of the g.c.d. and the 
substitution check imbedded in the program "Solve." The matrix of the homogeneous equations 
coefficients is read into the computer from magnetic tape no. 1. Other data to the program is 
provided by three data cards. The first card contains three numbers: the number of rows and 
columns of the matrix, and the number of primes supplied. The second and third cards contain 
the supplied primes. 

Functions used in this program are: 

"Inv" — computes the inverse of an integer t mod p which is required for solving the system 
mod p. 

"Irem" — IREM(^, t 2 , p) outputs an integer t satisfying t = £i£ 2 (mod p) , — p < t < p. Since 
the product t\t 2 is normally a double precision number, this routine is written in machine language. 

"Shift" — shifts a multiple-precision number to the right k positions until a binary 1 appears 
at the least significant position. On its left, k binary zeros are shifted in. Because of its multiple- 
precision character, this routine is also written in machine language. In the program listed below 
it appears directly after "Irem," before "Setpreci." 

"Mod" — MOD (t, p) transforms the single-precision integer t obtained as a result of an addition 
or subtraction of two numbers xu — (p — l)/2 ^ Xi ^ (p — l)/2, i = 1, 2, or of an IREM operation, 
to an integer £' satisfying £' = t(modp) and— (p — l)/2,^ t' ^ (p — l)/2. 

"Setpreci" — A Fortran subroutine which controls the multi-precision package which is written 
in machine language. Because this package is in a binary deck form, it is not included in the 
program listed below. Functions used from this package are: "If(Itest(£))" — for checking if a 
multiple precision integer t is negative zero or positive; "Multout" — for printing a multiple precision 
integer. 

PROGRAM SOLVE 

COMMON M1,M2,NZ,M,N,KP,KQ,KI,MA(111,120),IZ(120),NEW(9 ,9), 

1 NT(4,120),KT(20),IY(120)JSOL(120),KSOL(120)JS(112,9,6),KG(8) 

2 ,KE(8),KF(8),KO(8),NA(6),NB(6),NC(6),ND(6) 
TYPE MULTIPL6(6) A,B,C,KG,D,GCD,KC ,KSOL 
COMMON/SHIFTS/KK 

EQUIVALENCE (A,NA),(B,NB),(C,NC),(D,ND) 

READ55,M,N,NZ 

READ797,(KT(I),I=l,nz) 

REWIND 1 $ DO 987 I = 1 ,M 
987 READ TAPE 1,(MA(I,J),J= 1,N) 

PRINT 7,M,N,NZ $ DO 63 1= 1,M 
63 I Y(I) = I $ KR = $ PRINT 2 ,(KT(I),I = 1 ,NZ) $ DO 94 J = 1 ,N 
94 IZ(J) = J $ DO 200 L - 1 ,NZ $ KP = KT(L) $ KQ = KP/2. $ KD = MK = 1 

DO 75 1= l,M$DO 75 J= 1,N$IS = MA(I,J)-(MA(I,J)/KP)*KP$ IF(IS)76,75,75 
76 IS = KP + IS 
75 MA(I,J) = IS 

23 DO 122 K= MK,KR $ IF(IT.EQ.1)10,43 
8 PRINT 88,KP $ MS = 1 $ GO TO 201 
43 I1 = IY(K)$J1 = IZ(K)$KK=MA(I1J1) 
10 MK=MK+1$IF(KK- 1)8,33,34 

34 KD = IREM(KD,KK,KP) $ IM = INV(KK,KP) $ DO 32 JL = MK,N $ J2 = IZ(JL) 
32 MA(I1J2)=IREM(MA(I1J2),IM,KP) 
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33 DO 122 I = MK,M f 12 = IY(I) $ IL = MA(I2 Jl) $ IF(IL.EQ.0)122,31 

31 DOIJ = MK,N$J2 = IZ(J)$IS = MA(I2 J2) - IREM(MA(I1 J2),IL,KP)$IF(IS)72,1,1 

72 IS = KP + IS 

1 MA(I2J2) = IS 
122 CONTINUE 

D06I = MK,M$I1 = IY(I)$D06J = MK,N$J1 = IZ(J)$KK= MA(I1J1)$IF(KK.EQ.0)6,5 

5 KR = MK$IT = 1$IY(I) = IY(KR)$IY(KR) = I1$IZ(J) = IZ(KR)$IZ(KR) = J 1IGOT023 

6 CONTINUE $ IW = IT- 1 $ KF(IT)= KP $ IS= 1 $ KE(1)= 1 $ DO 70 1= 1,IW 
II = I + 1 $ IS = IREM(IS,KF(I),KP) 

70 KE(I1)=IS 

PRINT 51,KP,KD,IT $ IF(IT.EQ.1)47,29 
29 IS = KE(IT) $ IV = INV(IS,KP) $ NU = $ DO 944 KL = 1 ,IW 

KS = IREM(KO(KL),KE(KL),KP) $ IF(KS)971 ,979,979 
971 KS = KP + KS 

979 NU = NU + KS $ IF(NU - KP)944,974,974 
974 NU = NU-KP 

944 CONTINUE $ IS = KD - NU $ MS = KO(IT)= MOD(IREM(IS,IV,KP)) $ GO TO 95 
47 MS = 1 $ KO(l) = MOD(KD) $ KI = N - KR $ PRINT 56,KR,KI,(IY(I),I = 1 ,M) 

PRINT 58,(IZ(J),J = 1,N) 
95 DO26K=l,KI$DO26I=l,KR$Il=I-l$Kl=KR-Il$MX=0$DO38J=l,Il$K2=MK-J 
MX = MX-IREM(JSOL(K2),MA(IY(Kl),IZ(K2)),KP)$IF(MX)78,38,38 

78 MX = MX + KP 

38 CONTINUEIIS = MX - IREM(KD,MA(IY(K1),IZ(KR + K)),KP)$IF(IS)37 27,27 
37 IS = IS + KP 

27 JSOL(K1)=IS$IF(IT.EQ.1)28,90 

28 JS(Kl,K,l) = MOD(IS)$GOT0 26 

90 NU=0$DO44KL=l,IW$KS=IREM(JS(Kl,K,KL),KE(KL),KP)$IF(KS)71,79,79 

71 KS = KP + KS 

79 NU = NU + KS $ IF(NU - KP)44,74,74 
74 NU = NU-KP 

44 CONTINUE $ JS(K1 ,K,IT) = $ IS = IS - NU $ IF(IS.EQ.0)26,45 

45 MS= 1 $ JS(Kl,K,IT)=MOD(IREM(IS,IV,KP)) 
26 CONTINUE $ IT = IT +1 

201 REWIND 1 $ DO 202 I = 1 ,M $ RE AD TAPE 1,(MA(I,J)J = 1,N) 

202 CONTINUE $ IF(MS.EQ. 0)203 ,200 
200 CONTINUE $ GO TO 80 

203 PRINT 153 
KG(l)=l$DO710I = 2,IW 

710 KG(I) = KG(I-1)*KF(I-1) 

KC = KO(l) $ DO 720 1 = 2,IW 
720 KC = KC + KO(I)*KG(I) 

C = KC $ IV = 7 - 1 W $ IF(ITEST(C))751 ,80,752 

751 C = -C 

752 CALL SHIFT(NC(IV),IW) $ I5=KK 

DO 49 J = 1 JCI $ NOS = 15 $ A = C $ PRINT 100 J 

DO 50 1 = 1 ,KI $ KSOL(KR + 1) = 
50 CONTINUE $KSOL(KR + J) =KC 

DO40I = l,KR$B = JS(IJ,l)$JSOL(I)=l$DO42K = 2,IW 
42 B = B + JS(IJ,K)*KG(K) 

CALL MULTOUT(B,51) 

KSOL(I)=B 

IF(ITEST(B))715,716,717 
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716 JSOL(I) = 0$GOTO40 
715 B = -B$JS0L(I) = -1 

717 D0 718K=1,IW$I1 = 7-K 

718 JS(IJ,K)= NB(I1) $ CALL SHIFT(NB(IV),IW) $ IF(KK.GE.NOS)721,722 
722 NOS = KK 

721 D = A-B 

IF(ITEST(D))725,40,724 
725 D = -D 

724 CALL SHIFT(ND(IV),IW) $ A= B $ B = D $ GO TO 721 
40 CONTINUE $ GCD = A*2**NOS 

PRINT 669 

CALL MULTOUT(GCD,51) 

DO 730 1= 1,KR $ IF(JSOL(I).EQ.0)730,731 

731 D0 732K=1,IW$I1 = 7-K 

732 ND(I1) = JS(IJ,K) $ NV = D = D/GCD $ J2 = IZ(I) $ IF(JSOL(I).EQ. 1)733 ,734 
734 NV = -NV 

733 PRINT 538,(NT(K,J2),K=1,M1) 
PRINT 539,NV 

730 CONTINUE $ J2 = IZ(KR + J) $ NV = KC/GCD 
PRINT 538,(NT(KJ2),K= 1,M1) 
PRINT 539 ,NV 
DO61I=l,M$A = 0$DO52JL=l,N 

52 A = A + KSOL(JL)*MA(I,IZ(JL))$ IF(ITEST(A))82,61,82 
82 PRINT 53 J,I 

61 CONTINUE 
49 CONTINUE 

2 FORMAT(//50X,8H PRIMES.,//(I14,2X)) 

7 F0RMAT(1H1,3H M = ,I3,3H N = ,I3,4H NZ = ,I3,//,20X,7H MATRIX,//) 
51 FORMAT(1H1,50X6HPRIME = ,I14/I50,12H = DETERMINANT/150,1 1H 

ITERATIONS) 

53 FORMAT(1H1,20X,8H NOT YET,20X,2I20) 

55 FORMAT(23I5) 

56 FORMAT(//20X,11HFIRSTPRIME,5X,5HRANK=,I3,5X,8HNULLITY=,I3//40X,14 
1HORDER OF LINES, //(23Iff)) 

58 FORMAT(//,40X,16HORDER OF COLUMNS,//(23I5)) 

88 FORMAT(50XI14,23H THIS PRIME WAS DROPPED) 
100 FORMAT(1H1,//50X,I3///) 
153 FORMAT(/,20X,32HCONVERGENCE OBTAINED AND CHECKED,//) 

538 FORMAT(/,20X,10I3) 

539 FORMAT(I20, 8HX1 X2 X3,//) 
669 FORMAT(1H1,//50X,6H G.C.D//) 
797 FORMAT(5(I14,2X)) 

80 RETURN $ END 

FUNCTION INV(KX,KP) 

Kl =KP f K2 = KX $ Ml =0 $ M2 = 1 $ IF(K2 -1)2,4,5 

4 INV=K2$ RETURN 

6 M1=M2$M2 = M3$K1=K2$K2 = K3 

5 IQ = K1/K2 $ M3 = M1-M2*IQ $ K3 = K1-K2*IQ $ IF(K3 -1)2,7,6 

7 K2 = M3-(M3/KP)*KP $ IF(K2)8,2,4 

8 K2 = K2 + KP$GOT0 4 



2 INV = 
END 
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FUNCTION MOD(IX) 






COMMON M1,M2;M3,M4,M5,KP,KQ 




IF(IX)1,2,3 








2 MOD = IX#RETURI\ 


[ 




1 IF(KQ + IX)4,2,2 








4 MOD = KP + IX $ RETURN 




3 IF(KQ-IX)5,2,2 








5 MOD = IX-KP 








END 








IDENT 




IREM 




ENTRY 




IREM 


IREM 


SLJ 




** 




SIU 


1 


IR3 




LIU 


1 


*-l 




SIL 


1 


IR1 




INI 


1 


1 




SIL 


1 


IR2 




INI 


1 


1 


IR1 


SIL 


1 


IR3 




LIU 


1 


** 




LDA 


1 






LIL 


7 


IR1 


IR2 


MUI 


1 






LIU 


1 


** 




DVI 


1 






LLS 




48 


IR3 


ENI 


1 


** 




SLJ 




** 




END 








IDENT 




SHIFT 


' 


ENTRY 




SHIFT 


SHIFTS 


BLOCK 








COMMON 


NSHIFT 


SHIFT 


SLJ 




** 




SIU 


1 


EX 




SIL 


2 


EX 




LIU 


1 


SHIFT 




ENA 








STA 




= SNZ 




LDA 


1 






SAL 




N 




ARS 




24 




INI 


1 


1 


N 


SIU 


1 


EX + 2 




LIL 


1 


** 




INI 


1 


-1 




SAL 




BEGIN 




SAL 




CHANGE 




SAL 




LOOP 




INA 




1 




SAU 




STORE 


BEGIN 


ENI 


2 
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LDA 


1 


** 




AJP 




ZERO 




ENQ 






SH 


LRS 




1 




QJP 


M 


NEG 




INI 


2 


1 




SLJ 




SH 


NEC 


LLS 




1 




ENQ 


2 






STQ 




NSHIFT 




QJP 




EX 




UP 


1 


LOOP 


EXIT 


STA 


7 


BEGIN 


EX 


ENI 


1 


** 




ENI 


2 


** 




LDA 




NZ 




RAD 




NSHIFT 




SLJ 




** 


LOOP 


STA 




= ST 




LDA 


1 


** 




ENQ 








LRS 


2 






STA 




= ST1 




QRS 




1 




LDA 




T 




ADL 




= 03777777777777777 


STORE 


STA 


1 


** 




LDA 




Tl 




UP 


1 


LOOP 




SLJ 




EXIT 


ZERO 


SIL 


1 


= SNN ( 


ZEROl 


INI 


2 


47 




UP 


1 


CONT 




ENA 




-1 




STA 




NSHIFT 




SLJ 




EX 


CONT 


LDA 


7 


BEGIN 




AJP 




ZEROl 




SIL 


2 


= SNZ 




LIL 


2 


NN 


CHANGE 


LDA 


7 


BEGIN 




STA 


2 


** 




INI 


2 


-1 




UP 


1 


CHANGE 




ENA 






+ 


STA 


7 


CHANGE 




UP 


2 


* 




LIL 


1 


NN 



SLJ BEGIN 

END 
SUBROUTINE SETPRECI(NUM,LOG) 
COMMON / MULTINCM / N,IA(50),IC(50),ID(50) 
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' 



DATA (N=6) 
RETURN 
ENTRY ERR2 
PRINT 100,NUM,LOG 

100 FORMAT (46H ERROR IN ADDITION TYPE OVER-FLOW CALL FROM 
1,016, ,18H FIRST LOCATION ,016) 

STOP 

ENTRY ERR1 
PRINT 101,NUM,LOG 

101 FORMAT (49H ERROR IN DIVISION TYPE ZERO-DIVISOR CALL 
1FROM ,016, ,16H FIRST LOCATION ,016 ) 

STOP 

ENTRY PRINTOUT 

WRITE (NUM.102) (ID(K),K=1,N) 

102 FORMAT (5X,5(I14,X)) 
END 

DATA 

111 120 10 
7908189600581 7908189600583 7908189600587 7908189600589 7908189600593 

7908189600599 7908189600601 7908189600607 7908189600611 7908189600613 
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